/**
  EG4Sim project, Anyang Normal University && IMP-CAS
  \class EG4DetConstruct.h
  \brief To construct detectors for G4 sim
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2019-01-21
  \date 2023-11-07 last modified
  \copyright Copyright (c) 2023-2024 IMP-CAS with LGPLv3 LICENSE
*/

#include "G4SDManager.hh"
#include "G4NistManager.hh"
#include "G4Material.hh"
#include "G4FieldManager.hh"
#include "G4LogicalVolume.hh"
#include "G4TransportationManager.hh"
#include "G4PVPlacement.hh" // inherited from G4VPhysicalVolume
#include "G4SystemOfUnits.hh" // units imported from CLHEP library
#include "G4UserLimits.hh"
#include "G4VisAttributes.hh"
#include "G4Colour.hh"
#include "G4Box.hh"

#include "EG4DetConstruct.h"
#include "EG4ParaManager.h"
#include "EG4DC.h"
#include "EG4DC6.h"
#include "EG4TA.h"
#include "EG4TA0.h"
#include "EG4DCTaD.h"
#include "EG4Dipo.h"
#include "EG4NW.h"
#include "EG4T0.h"
#include "EG4T1.h"
#include "EG4TOFW.h"
#include "EG4TOFW30.h"
#include "EG4MagneticField.h"
#include "EG4DCSD.h" // sensitive detectors for drift chambers
#include "EG4PlaSD.h" // sensitive detectors of plastic scintillators
#include "EG4NWSD.h"

using std::cout;
using std::endl;

// #define DCArrM // define the auxiliary dc array

EG4DetConstruct::EG4DetConstruct() : G4VUserDetectorConstruction(), fT0(),
fDCm(0), fDC6m(0), fTOFWm(0){}
EG4DetConstruct::~EG4DetConstruct(){}

void EG4DetConstruct::readpos(const Node &f, double *p, double *a){
  for(int i = 0; i < 3; i++){
    p[i] = f["pos"][i].as<double>();
    a[i] = f["ang"][i].as<double>();
  } // end for over i
} // end member function readpos

G4VPhysicalVolume *EG4DetConstruct::Construct(){
  fpm = EG4ParaManager::Instance();
  const Node &dp = fpm->DetPosPara();
  char name[128]{};

  ConstructMaterial(); // construct all the material needed globally
  // the world
  // fWorldLV = new G4LogicalVolume(new G4Box("World", 3.*m, 1.*m, 15.*m), fVacuum, "World");
  fWorldLV = new G4LogicalVolume(new G4Box("World", 3.5*m, 1.*m, 15.*m), fAir, "World");
  fWorld = new G4PVPlacement(0, G4ThreeVector(), fWorldLV,
    "World", nullptr, false, 0, true);
  G4VisAttributes *vat = new G4VisAttributes(false, G4Colour(1., 1., 1.)); // white
  fWorldLV->SetVisAttributes(vat);

  // T0, i.e., Tstart
  // fT0 = new EG4T0();
  // new G4PVPlacement(nullptr, G4ThreeVector(0., 0., -3000.),
  //   fT0->LV(), "T0", fWorldLV, false, 0, true);

  // T1, i.e., Tstop
  double pT1[3]{}, aT1[3]{};
  readpos(dp["tstop"], pT1, aT1);
  fT1 = new EG4T1();
  new G4PVPlacement(nullptr, G4ThreeVector(pT1[0], pT1[1], pT1[2]),
    fT1->LV(), "tstop", fWorldLV, false, 0, true);

  // the primary target
  // fTA0 = new EG4TA0();
  // new G4PVPlacement(nullptr, G4ThreeVector(0., 0., -3100.*mm), fTA0->LV(),
  //   "TA", fWorldLV, false, 0, true);

  // the reaction target
  // fTA = new EG4TA();
  // z of 32Q01: (the plane of the yellow magnet iron wall): -3205.82
  // new G4PVPlacement(nullptr, G4ThreeVector(0., 0., -3205.82*mm), fTA->LV(), // 3205.82
  //   "TA", fWorldLV, false, 0, true);
  // new G4PVPlacement(nullptr, G4ThreeVector(0., 0., -1806.39*mm), fTA->LV(),
  //   "TA", fWorldLV, false, 0, true);

  // DCTaD
  fDCTaD = new EG4DCTaD();
  // position measured @Dec.12, 2023
  // position measured @Dec.12, 2023 - z of the X0 wire layer
  double pDCTaD[2][3]{}; // center of the X0 wire layer
  double aDCTaD[2][3]{}; // in degree
  const G4ThreeVector x0DCTaD(0., 0., -17.5); // center w.r.t. x0
  for(int i = 0; i < 2; i++){
    sprintf(name, "dcTaD%d", i);
    readpos(dp[name], pDCTaD[i], aDCTaD[i]);
    auto dcRot = new G4RotationMatrix(); dcRot->rotateZ(-aDCTaD[i][2]*deg);
    dcRot->rotateX(-aDCTaD[i][1]*deg); dcRot->rotateY(-aDCTaD[i][0]*deg);
    // x, z: the geometrical center of the DC box
    auto pos = (*dcRot)*x0DCTaD + G4ThreeVector(pDCTaD[i][0], pDCTaD[i][1], pDCTaD[i][2]);
    new G4PVPlacement(dcRot, pos, fDCTaD->LV(), name, fWorldLV, false, i, true);
    cout << name << pos << endl; // DEBUG
  } // end for over i

  // the dipo magnet //
  fDipo = new EG4Dipo();
  fDipo->SetColor(0.9, 0.9, 0.9); // light gray
  new G4PVPlacement(0, G4ThreeVector(), fDipo->LV(), "Dipo", fWorldLV, false, 0, true);
  fDipo->LV()->SetVisAttributes(vat);

  // DC
  fDC = new EG4DC();
  // position measured @Dec.12, 2023 - z of the X0 wire layer
  double pDC[3][3]{}; // center of the X0 wire layer
  double aDC[3][3]{}; // in degree
  const G4ThreeVector x0DC(0., 0., 5.); // center w.r.t. x0
  for(int i = 0; i < 3; i++){
    sprintf(name, "dc%d", i);
    readpos(dp[name], pDC[i], aDC[i]);
    auto dcRot = new G4RotationMatrix(); dcRot->rotateZ(-aDC[i][2]*deg);
    dcRot->rotateX(-aDC[i][1]*deg); dcRot->rotateY(-aDC[i][0]*deg);
    // x, z: the geometrical center of the DC box
    auto pos = (*dcRot)*x0DC + G4ThreeVector(pDC[i][0], pDC[i][1], pDC[i][2]);
    new G4PVPlacement(dcRot, pos, fDC->LV(), name, fWorldLV, false, i, true);
    cout << name << pos << endl; // DEBUG
  } // end for over i

  // TOFW
  fTOFW = new EG4TOFW();
  double pTOFW[3]{}, aTOFW[3]{};
  readpos(dp["tofwB"], pTOFW, aTOFW);
  auto tofwRot = new G4RotationMatrix(); tofwRot->rotateZ(-aTOFW[2]*deg);
  tofwRot->rotateX(-aTOFW[1]*deg); tofwRot->rotateY(-aTOFW[0]*deg);
  new G4PVPlacement(tofwRot, G4ThreeVector(pTOFW[0], pTOFW[1], pTOFW[2]),
    fTOFW->LV(), "tofwB", fWorldLV, false, 0, true);
  fTOFW->LV()->SetVisAttributes(vat);

  // the Newtron Wall (NW)
  fNW = new EG4NW();
  new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 13066.),
    fNW->LV(), "NW", fWorldLV, false, 0, true);
  fNW->LV()->SetVisAttributes(vat);
  fNW->VetoLV()->SetVisAttributes(vat);
  fNW->SciLayerLV()->SetVisAttributes(vat);
  fNW->FeLayerLV()->SetVisAttributes(vat);
  fNW->FeStripLV()->SetVisAttributes(vat);

  //-------- THE OTHER DC ARM ---------//
  // DCm  
  double pDCm[3][3]{}; // center of the X0 wire layer
  double aDCm[3][3]{}; // in degree
  for(int i = 0; i < 3; i++){
    sprintf(name, "dcm%d", i);
    if(dp[name]){
      if(!fDCm) fDCm = new EG4DC();
      if(!fDC6m) fDC6m = new EG4DC6();
      readpos(dp[name], pDCm[i], aDCm[i]);
      auto dcRot = new G4RotationMatrix(); dcRot->rotateZ(-aDCm[i][2]*deg);
      dcRot->rotateX(-aDCm[i][1]*deg); dcRot->rotateY(-aDCm[i][0]*deg);
      // x, z: the geometrical center of the DC box
      auto pos = (*dcRot)*x0DC + G4ThreeVector(pDCm[i][0], pDCm[i][1], pDCm[i][2]);
      new G4PVPlacement(dcRot, pos, 2 == i ? fDC6m->LV() : fDCm->LV(), name,
        fWorldLV, false, i, true);
      cout << name << pos << endl; // DEBUG
    } // end if
  } // end for over i

  // TOFWm
  double pTOFWm[3]{}, aTOFWm[3]{};
  if(dp["tofwm"]){
    if(!fTOFWm) fTOFWm = new EG4TOFW30();
    readpos(dp["tofwm"], pTOFWm, aTOFWm);
    auto tofwRotm = new G4RotationMatrix(); tofwRotm->rotateZ(-aTOFWm[2]*deg);
    tofwRotm->rotateX(-aTOFWm[1]*deg); tofwRotm->rotateY(-aTOFWm[0]*deg);
    new G4PVPlacement(tofwRotm, G4ThreeVector(pTOFWm[0], pTOFWm[1], pTOFWm[2]),
      fTOFWm->LV(), "tofwB", fWorldLV, false, 0, true);
    fTOFWm->LV()->SetVisAttributes(vat);
  } // end if

  // --- the tool room --- //
  // auto trlv = new G4LogicalVolume(new G4Box("toolroom", 1.*m, 1.*m, 1.*m),
  //   fAir, "toolroom");
  // auto trrot = new G4RotationMatrix(); trrot->rotateY(30.*deg);
  // new G4PVPlacement(trrot, G4ThreeVector(4904.18*mm, 0., 2439.18*mm),
  //   trlv, "toolroom", fWorldLV, false, 0, true);

  return fWorld;
} // end of member function Construct

inline void setSD(G4SDManager *m, G4LogicalVolume *l, G4VSensitiveDetector *s){
  m->AddNewDetector(s);
  l->SetSensitiveDetector(s);
}
void EG4DetConstruct::ConstructSDandField(){
  // sensitive detectors //
  G4SDManager *sdm = G4SDManager::GetSDMpointer();
  if(fT0) setSD(sdm, fT0->LV(), new EG4PlaSD("T0_SD"));
  setSD(sdm, fT1->LV(), new EG4PlaSD("T1_SD"));
  setSD(sdm, fDCTaD->WireLV(), new EG4DCSD("DCTaD_SD"));
  setSD(sdm, fDC->WireLV(), new EG4DCSD("DC_SD"));
  setSD(sdm, fTOFW->StripLV(), new EG4PlaSD("TOFWStrip_SD"));

  if(fDCm) setSD(sdm, fDCm->WireLV(), new EG4DCSD("DCm_SD"));
  if(fDCm) setSD(sdm, fDC6m->WireLV(), new EG4DCSD("DC6m_SD"));
  if(fDCm) setSD(sdm, fTOFWm->StripLV(), new EG4PlaSD("TOFWmStrip_SD"));

  setSD(sdm, fNW->SciStripLV(), new EG4NWSD("NWSciStrip_SD"));
  setSD(sdm, fNW->Sci10LV(), new EG4NWSD("NWFeStrip_SD"));
  setSD(sdm, fNW->VetoStripLV(), new EG4PlaSD("NWVetoStrip_SD"));

  // magnetic field //
  EG4MagneticField *mag = new EG4MagneticField();
  auto mf = G4TransportationManager::GetTransportationManager()->GetFieldManager();
  mf->SetDetectorField(mag);
  // create the objects which calculate the trajectory for a pure magnetic field
  mf->CreateChordFinder(mag);
  fDipo->LV()->SetFieldManager(mf, true); // true: force to all daughters
} // end member function CreateSDandField

// construct needed material for global use
void EG4DetConstruct::ConstructMaterial(){
  // get Lanzhou local air
  static const double t = 293.15*kelvin, p = 0.841*atmosphere; // local for Lanzhou
  G4NistManager *nist = G4NistManager::Instance();
  fAir = nist->ConstructNewGasMaterial("airLanzhou", "G4_AIR", t, p);
  fVacuum = new G4Material("Galactic", 1., 1.01*g/mole,
    CLHEP::universe_mean_density, kStateGas, 0.1*kelvin, 1.e-19*pascal);

  // dc gas mxiture //
  fAr = nist->ConstructNewGasMaterial("arLanzhou", "G4_Ar", t, p);
  fCO2 = nist->ConstructNewGasMaterial("co2Lanzhou", "G4_CARBON_DIOXIDE", t, p);
  const double arVolF = 0.8, co2VolF = 0.2;
  const double dcDen = fAr->GetDensity()*arVolF + fCO2->GetDensity()*co2VolF;
  fDCGas = new G4Material("c20Lanzhou", dcDen, 2, kStateGas, t, p);
  fDCGas->AddMaterial(fAr, arVolF*39.95/(arVolF*39.95+co2VolF*44.01));
  fDCGas->AddMaterial(fCO2, co2VolF*44.01/(arVolF*39.95+co2VolF*44.01));

  // use 5052 Aluminum alloy for the timing being
  fAl5052 = new G4Material("Al5052", 2.68*g/cm3, 3);
  G4Element *Al = nist->FindOrBuildElement("Al");
  G4Element *Mg = nist->FindOrBuildElement("Mg");
  G4Element *Cr = nist->FindOrBuildElement("Cr");
  fAl5052->AddElement(Al, 97.266*perCent);
  fAl5052->AddElement(Mg, 2.252*perCent);
  fAl5052->AddElement(Cr, 0.482*perCent);

  //----- the epoxy resin plates widely used in PCB -----//
  G4Element *Si = nist->FindOrBuildElement("Si");
  G4Element *O = nist->FindOrBuildElement("O");
  G4Element *H = nist->FindOrBuildElement("H");
  G4Element *C = nist->FindOrBuildElement("C");
  // quartz
  G4Material *SiO2 = new G4Material("quartz", 2.200*g/cm3, 2);
  SiO2->AddElement(Si, 1);
  SiO2->AddElement(O , 2);
  // epoxy (for FR4)
  G4Material* epoxy = new G4Material("epoxy", 1.2*g/cm3, 2);
  epoxy->AddElement(H, 2);
  epoxy->AddElement(C, 2);
  // FR4 (Glass + epoxy)
  fFR4 = new G4Material("FR4", 1.86*g/cm3, 2);
  fFR4->AddMaterial(SiO2, 52.8*perCent);
  fFR4->AddMaterial(epoxy, 47.2*perCent);
} // end member function ConstructMaterial
